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Abstract 

We  consider  goal-oriented  a  posteriori  error  estimators  for  the  evaluation  of  the  er¬ 
rors  on  quantities  of  interest  associated  with  the  solution  of  geometrically  nonlinear 
curved  elastic  rods.  For  the  numerical  solution  of  these  nonlinear  one-dimensional 
problems,  we  adopt  a  B-spline  based  Galerkin  method,  a  particular  case  of  the 
more  general  Isogeometric  Analysis.  We  propose  error  estimators  using  higher  or¬ 
der  “enhanced”  solutions,  which  are  based  on  the  concept  of  enrichment  of  the 
original  B-spline  basis  by  means  of  the  “pure”  fc-refinement  procedure  typical 
of  Isogeometric  Analysis.  We  provide  several  numerical  examples  for  linear  and 
nonlinear  output  functionals,  corresponding  to  the  rotation,  displacements  and 
strain  energy  of  the  rod,  and  we  compare  the  effectiveness  of  the  proposed  error 
estimators. 

Keywords:  Geometrically  nonlinear  rods;  Isogeometric  Analysis;  B-spline  basis;  goal- 
oriented  a  posteriori  error  estimation;  error  estimator. 


1  Introduction 

The  study  of  large  deflections  of  thin  beams,  or  rods,  has  received  a  growing  interest 
in  many  engineering  and  science  problems.  Examples  of  these  problems  include  framed 
structures  [36,  37],  compliant  mechanisms  [11,  44]  and  nanoscale  structures  [45]. 

The  simplest  geometrically  nonlinear  bending  theory  of  planar  elastic  rods  is  the  Elas- 
tica  theory ,  mathematically  formulated  by  L.  Euler  in  1744  [31].  In  this  theory  a  rod 
is  thought  of  as  an  inextensible  line  of  particles  which  resists  bending  according  to  a 
law  given  by  a  linear  constitutive  relation;  no  restrictions  on  the  magnitude  of  dis¬ 
placements  or  angles  of  rotation  are  considered.  Since  then,  several  variants  have  been 
proposed,  in  particular,  theories  including  dynamical  effects,  extensibility  of  the  rod, 
shear  deformation,  plasticity,  follower  loads,  etc.,  see  e.g.  [3]  and  the  references  therein. 
Also,  different  approaches  have  been  considered  in  the  literature  for  the  analysis  of 
elastic  thin  beams:  (i)  the  elliptic  integral  approach  first  proposed  by  [9],  which  gives 
closed-form  solutions  only  for  simple  loading  cases  and  boundary  conditions  [33] ,  (ii)  the 
numerical  integration  approach  with  iterative  shooting  techniques  (e.g.  [29,  32,  13]),  and 
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(iii)  the  incremental  Finite  Element  method  with  Newton-Raphson  iteration  techniques 
[43,  22,  21], 

Of  these  approaches,  the  Finite  Element  method  ( e.g .  [14,  40,  25])  is  indeed  the  most 
popular  approach,  mainly  due  to  its  versatility  to  the  analysis  of  problems  with  complex 
topologies  and  geometries.  As  a  result,  numerous  geometrically  nonlinear  planar  beam 
elements  have  been  developed  over  the  past  few  decades.  Among  others,  see  e.g.  [17, 
19,  18,  24,  41,  27,  36,  47,  12,  30,  23,  37,  42], 

However,  the  Finite  Element  method  only  provides  approximate  solutions  whose 
quality  depends  on  the  discretization  procedure.  The  problem  of  how  to  measure  the 
quality  of  these  numerical  solutions,  assessing  the  accuracy  of  the  Finite  Element  ap¬ 
proximations,  is  essential  for  a  reliable  application  of  these  problems. 

Nevertheless,  if,  on  the  one  hand,  error  estimation  offers  no  major  challenges  in  the 
linear  analysis  of  one-dimensional  beam  problems,  since  it  is  possible  to  employ  certain 
elements  which  can  model  the  linear  problem  exactly,  even  with  one  element  per  member, 
on  the  other  hand,  approaches  based  on  a  posteriori  error  estimation  [2]  for  geometrically 
nonlinear  beams  have  received  very  little  attention.  Exceptions  to  this  are  the  works 
presented  in  [28]  and  [34]. 

Furthermore,  although  the  error  estimation  on  the  solution  of  a  problem  represents 
an  important  aspect  of  the  numerical  approximation,  it  is  also  important  to  properly 
evaluate  quantities  of  interest  associated  with  such  solution,  often  referred  to  as  output 
functionals  (see  e.g.  [4]).  For  the  rod  problem,  these  quantities  of  interest,  can  be 
defined  by  local  quantities,  such  as,  for  instance,  displacements,  rotations,  forces  and 
bending  moments,  or  by  global  quantities,  such  as  the  total  strain  energy.  In  this 
context,  a  posteriori  error  estimation  approaches,  are  often  referred  to  as  goal-oriented 
approaches,  as  they  are  particularly  suited  for  estimating  the  error  on  quantities  of 
interest,  rather  than  the  solution  itself;  the  introduction  of  the  dual  (adjoint)  problem 
allows  the  quantification  of  the  sensitivity  of  the  output  functionals  with  respect  to  the 
perturbations  on  the  solution.  Several  contributions  have  been  made  in  this  field,  with 
special  emphasis  on  the  Finite  Element  method,  for  both  linear  and  nonlinear  problems, 
see  e.g.  [2,  8,  35,  20,  5]. 

It  is  the  main  goal  of  this  work  to  propose  and  discuss  goal-oriented  a  posteriori 
error  estimators  for  the  evaluation  of  the  errors  on  quantities  of  interest  involving  the 
solution  of  geometrically  nonlinear  rods  problems.  With  this  aim,  we  consider  a  numer¬ 
ical  approximation  scheme  based  on  a  B-spline  basis  [38]  for  the  approximation  space. 
This  choice  represents  a  particular  case  of  the  Isogeometric  Analysis  method  [26,  15],  a 
Galerkin  approximation  method  based  on  the  isoparametric  concept  in  which  the  basis 
of  the  approximation  space  is  the  same  as  the  one  used  to  represent  the  geometry  and 
which  can  be  represented  by  B-splines,  NURBS  (Non-Uniform  Rational  B-splines)  [38] 
or  eventually  T-splines  [6].  High  order,  Ck  globally  continuous  basis  can  be  generated 
while  maintaining  the  exactness  of  the  representation  of  the  geometry  in  the  analysis; 
moreover,  the  basis  possesses  useful  properties  for  the  approximation  of  solutions  inde¬ 
pendently  of  the  geometrical  construction.  Error  analysis  methods  [7]  and  a  posteriori 
error  estimates  with  T-splines  [16]  have  already  been  proposed  in  the  literature,  also  in 
the  goal-oriented  framework  [46] . 

The  proposed  error  estimators  are  constructed  on  a  B-spline  basis  in  which  the  “en¬ 
hanced”  approximated  solutions,  typically  employed  to  evaluate  the  estimators,  are  not 
only  generated  by  means  of  a  mesh  refinement  strategy,  as  proposed  in  [46]  for  two- 
dimensional  B-splines  and  T-splines,  but  also  by  means  of  higher  order  approximations. 
In  particular,  we  build  the  “enhanced”  higher  order  B-spline  basis  by  performing  one 
step  of  the  “pure”  fc-rehnement  procedure  typical  of  Isogeometric  Analysis  [15]  (order 
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Figure  1:  Cantilever  rod  model. 


elevation  without  internal  knots  insertion);  in  this  manner,  an  higher  order  B-spline  ba¬ 
sis  with  increased  global  continuity  is  obtained  by  introducing  only  an  additional  degree 
of  freedom  with  respect  to  the  original  one.  As  a  result,  the  corresponding  estimators 
and,  therefore,  the  errors  associated  with  the  outputs,  can  be  evaluated  at  a  relatively 
small  computational  cost.  The  case  of  the  classical  Finite  Element  method  with  the 
standard  Lagrangian  polynomial  linear  basis  is  obtained  and  discussed  as  a  particular 
case  of  the  B-spline  basis  of  order  one. 

As  the  proposed  error  estimators  are  formulated  in  a  general  framework,  they  can  be 
straightforwardly  applied  to  the  analysis  of  other  nonlinear  problems,  even  multidimen¬ 
sional. 

The  paper  is  organized  as  follows.  In  Sec. 2  we  introduce  the  boundary  value  problem; 
in  Sec. 3  we  reformulate  the  problem  in  an  abstract  setting  by  introducing  its  correspond¬ 
ing  weak  form  and  recall  goal-oriented  a  posteriori  error  estimates.  In  Sec. 4  we  describe 
the  B-spline  based  approximation  and  propose  the  error  estimators.  In  Sec. 5  we  pro¬ 
vide  numerical  tests  for  the  rod  problem  and  discuss  the  obtained  results.  Conclusions 
follow. 


2  Boundary— value  problem 

We  confine  attention  to  the  statics  of  cantilever  rods  having  a  planar  initially  curved 
reference  configuration  subjected  to  a  pair  of  concentrated  loads  Px  and  Py  and  a  bend¬ 
ing  moment  M  applied  at  the  tip,  as  represented  in  Fig.l.  The  governing  differential 
equation  of  these  type  of  rods  is  given  as  follows  [31]: 

El-j-^  +  Px  sin(0  +  6q)  +  Py  cos (0  +  9q)  =  0.  (1) 

where  s  £  fl  =  (0,  L)  represents  the  arch  length  of  the  rod  (curvilinear  abscissa)  in  the 
undeformed  configuration,  with  L  the  length  of  the  rod;  9q  =  9q(s)  stands  for  the  initial 
slope  angle  of  the  rod,  9  =  0(s)  represents  the  rotation  angle  of  the  normal  to  rod  axis, 
and  El  is  the  bending  stiffness  of  the  beam,  which  we  assume  as  constant,  being  E  the 
Young  modulus  and  I  the  cross-sectional  moment  of  inertia.  The  boundary  conditions 
for  the  differential  equation  (1)  are: 


0(0)  =  0, 


d0 

EITs 


=  M. 


s—L 


(2) 
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As  it  is  well  known,  the  boundary-value  problem  defined  by  Eqs.(l)  and  (2)  may  exhibit 
multiple  solutions,  very  often  referred  in  the  literature  to  as  elastica  shapes.  Each 
shape  corresponds  to  a  critical  point  of  the  total  potential  energy  functional  defined  by 

iy-)  :  V^K: 

np(0)  :=  17(0)  +  W(6)  (3) 

where  U  and  W  represent  the  strain  (internal)  and  the  external  energies  of  the  rod 
element,  respectively: 

ui-e)--=  <4> 


W(0)  := 


/  (Px  cos (0  +  0O )  —  Py  sin(0  +  0O))  ds 
J  n 

+  Px(x o(0)  -  x0(L)) 


Py(yo(0)-yo{L))-Me(L),  (5) 


with  V  the  kinematically  admissible  (functional)  space  defined  as: 

V  :=  {0  €  :  0(0)  =0}, 


(6) 


where  H1  (fi)  represents  a  standard  Hilbert  space  [1];  the  pairs  (xo(0),  2/o (0))  and  (xq(L),  yo(L)) 
represent  the  initial  coordinates  of  the  cantilever  at  s  =  0  and  s  =  L,  respectively. 

The  horizontal  and  vertical  displacements  of  the  rod,  herein  denoted  by  u  =  u(s)  and 
v  =  v(s),  respectively  (see  Fig.l),  are  shown  to  obey  the  following  differential  equations: 


du 

ds 

dv 

ds 


cos  (0  +  0O ) 
sin(0  +  0O ) 


dx  o 
ds  ’ 
dyo 
ds 


where: 

dx  o 
ds 
dyo 
ds 

together  with  the  boundary  conditions: 


cos(0o), 

sin(0o), 


u(0)  =  n(0)  =  0. 


(7) 

(8) 

(9) 


We  are  interested  in  evaluating  quantities  of  interest  associated  with  the  solution  0(s) 
of  problem  (1)— (2) ,  in  particular,  the  strain  energy  17(0)  defined  by  (4),  the  rotation  at 
the  tip  of  the  cantilever  rods  0(7),  and  the  horizontal  and  vertical  displacements  at  the 
tip  of  the  cantilever  rods,  u(L)  and  v{L).  The  displacements  can  be  obtained  from  the 
integration  of  equations  (7)  and  (9)  as: 


u(L)  =  f  (cos(0  +  0O )  —  cos(0o))  ds. 
Jn 

v(L)  =  f  (sin(0  +  0O )  —  sin(0o))  ds. 
Jn 


(10) 


3  Goal— oriented  a  posteriori  error  estimation 

In  this  Section  we  recall  and  discuss  the  so-called  goal-oriented  analysis;  with  this  aim, 
we  reformulate  the  problem  discussed  in  the  previous  Section  in  an  abstract  setting  by 
introducing  its  weak  form  and  the  adopted  Galerkin  approximation  scheme. 
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3.1  Weak  form 

The  weak  form  associated  with  the  primal  boundary-value  problem  introduced  in  the 
preceding  Section  reads: 

find  9  G  V  :  a{9){4>)  =  V</>  €  V,  (11) 

with  the  form  a(-)(-)  :  V  x  V  — >  R.  semilinear  in  the  second  argument  (nonlinear  in  the 
first)  and  differentiable  in  Frechet  sense,  and  the  functional  /(•)  :  V  — >  R.  linear  and 
continuous.  We  assume  that  suitable  hypothesis  for  the  existence  and  local  uniqueness 
of  the  solutions  of  problem  (11)  hold.  Additionally,  in  view  of  the  error  estimation,  it 
is  convenient  to  introduce  the  primal  residual  Rpr(-)(-)  :  V  x  V  — >•  R  as: 

Rr{9)(<l>):=f{<l>)-a(0){<t>).  (12) 

Associated  with  the  solution  of  problem  (11),  we  desire  to  evaluate  quantities  of  interest 
qi  represented  by  output  functionals: 

qi  =  U(9)  i  =  (13) 

with  li(9)  :  V  — >•  R  suitable  differentiable  linear  or  nonlinear  functionals. 

By  introducing  the  finite  dimensional  approximation  space  Vh  C  V,  the  Galerkin  ap¬ 
proximated  primal  problem  reads: 

find  9h  €  Vh  ■  a(9h){<f>h)  =  f{<j)h)  V<t>h  &  Vh,  (14) 


with  the  associated  approximated  output  functionals  q^: 

qi,h  =  k(9h )  i  =  l,...N.  (15) 

If  we  refer  to  the  problem  (1) — (2),  we  have  that  V  =  {v  £  Hx(fl)  :  v(0)  =  0}  and: 
rd9  dip 


a{9){4>)  = 


El 


ds  ds 


—  Px  sin(0  +  6o)<j)  —  Pv  cos (9  +  6q)(/)  I  ds 


(16) 


m  =  M<KL). 

We  are  interested  in  evaluating  the  following  output  functionals 

qi  =  h(9):=U(9),  q2  =  l2(9)  :=  9(L ), 

<?3  =  h{9)  :=  u(L),  qi  =  U(9)  :=  v(L), 


(17) 


where  l\{9)  is  quadratic  in  9 ,  h(9)  is  linear  in  9  and  h(9)  and  li{9)  are  nonlinear 
trigonometric  functionals.  Note  that,  in  general,  <72  =  9{sq)  can  be  represented  as 

follows  I2 (6)  =  /  9(s)S(s  —  So)ds,  with  S(s)  the  delta  Dirac  function  and  So  6  D- 
Jn 

3.2  A  posteriori  error  estimation 

A  posteriori  error  estimation  in  the  goal-oriented  framework  represents  a  well  estab¬ 
lished  tool  for  errors  associated  with  output  functionals,  see  e.g.  [8,  35,  20,  5],  especially 
in  the  linear  case.  In  this  Section  we  recall  results  for  nonlinear  problems  in  view  of  the 
definition  of  the  error  estimators. 
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We  introduce  the  dual  variables  Zi  £  V,  representing  the  solutions  of  the  following 
dual  problems: 

find  Zi&V  :  a\e){z^)  =  l'i{e){4>)  V0  £  V,  i  =  l,...,N,  (18) 

which  depend  on  the  primal  solution  9  £  V  (11)  and  on  the  output  functionals  IfiO)  (13). 

The  notations  a'(-)(-,  •)  :  Vx  Vx  and  /((•)(■)  :  V  x  V  — X  R  indicate  the  first 

Frechet  differentials  in  9  of  the  form  a(9)(-)  and  functionals  U{9).  Note  that,  for  a  given 
9  £  V,  the  form  a'(9)(-,  •)  :  V  x  V  — >  K  and  the  functionals  l[{9){-)  :  V  — >  1  are  linear; 
it  follows  that  the  dual  problems  (18)  are  linear  with  respect  to  the  dual  variables. 
Similar  considerations  follow  for  higher  order  differentials.  Additionally,  it  is  convenient 
to  introduce  the  dual  residuals  Rfu(-)(-,  •)  :  V  x  V  x  V  — >  R  as: 

i  =  l,...,N.  (19) 

The  Galerkin  approximation  of  the  dual  problems  (18)  reads: 

find  zi:h  £  Vh  ■  a'{9h)(zh,(f>h)  =  l'l(9h)((j)h)  \/<f)h  £  Vh,  i  =  l,...,N,  (20) 

with  9h  £  Vh  the  solution  of  the  primal  problem  (14). 

The  goal  of  the  a  posteriori  error  estimation  consists  in  properly  estimate  the  errors 
on  the  output  functionals 

\Qi~Qi,h\  i  =  (21) 

once  the  approximated  primal  9h  £  V  is  provided.  With  this  aim,  we  recall  from  [8] 
the  following  Propositions,  in  which  the  Galerkin  orthogonality  property  has  been  taken 
into  account: 


Proposition  3.1  For  the  Galerkin  approximated  primal  (14)  and  dual  (20)  problems, 
we  have: 

qi  -  qi,h  =  Sqi+lli  i=l,...,N,  (22) 

with: 

£qi  =  £qi(0,Zi)  :=  l-Rpr{9h){Zi)  +  ^Rtu(9h)(Zi,h,9),  (23) 

the  remainder  terms: 


Tlr  := 


i  r 1 

-  /  {IT {Oh  +  tepr)(epr ,  epr ,  epr )  -  a"'{9h  +  tepr){epr ,  epr,  epr,  zh  +  tefu) 

^  Jo 


-3 a"(9h  +  tepr)(epr,  epr ,  ef )}  t(t  -  1  )dt 


(24) 


and  the  primal  and  dual  errors  defined  as  epr  :=  9  —  9h  and  edu  :=  Zi  —  Zith- 


Proposition  3.2  Let  us  assume  that  9h  £  Vh  C  V  is  an  improved  (enhanced)  approxi¬ 
mated  solution  with  respect  to  9h  £  Vh,  then: 


qi  qi,h  —  £qi  V  'R'i  i  1 ,  •  •  • ,  -N, 

with: 

£qi=£qi(zi)  ~Rpr(9h)(zi), 

and  the  reminder  terms  such  that: 


’R, 


<  max  | /''(^)(?r,  epr)  -  a"(^)(^,  epr ,  z)\ , 

,8h,8h 


(25) 

(26) 

(27) 


being  9 ,  Oh,  9h  the  “triangle”  spanned  in  V  by  the  functions  9,  Oh  and  Oh  and  epr  :=  0—0h- 
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Let  us  observe  that: 


1.  the  a  posteriori  error  estimates  (22)  and  (25)  are  not  computable,  since  they 
involve  the  exact  solutions  9  £  V  and  Zi  £  V; 

2.  the  remainder  terms  IZi  in  Eq.(24)  are  null  for  at  most  quadratic  functionals  h(9) 

and  bilinear  form  =  a(9,  •); 

3.  in  the  case  of  linear  functionals  U(9)  and  bilinear  form  a(0)(-)  =  a(9,  •),  the  re¬ 
mainder  terms  74  and  IZi  hr  Eqs.(24)  and  (27)  are  null  and  the  errors  qt  —  qi^h  = 

£qi  =  £qi  i 

4.  if  the  differentials  of  the  form  a(9)(-)  and  functionals  li{9)  can  be  bounded  for 
9  £  V,  the  remainder  terms  IZi  and  7 Zi  converge  to  zero  with  order  3  and  2, 
respectively,  in  the  errors  epr  and  ef“; 

5.  for  the  rod  problem  under  consideration  (see  Eqs.(16)  and  (17)),  the  remainder 
terms  IZi  and  IZi  can  be  bounded,  since  the  functionals  involved  are  at  most 
quadratic  or  depend  on  trigonometric  functions  which  are  C°°  continuous  and 
bounded;  similar  properties  hold  for  the  form  a(9)(-); 

6.  the  terms  £qi  and  £qi  in  the  estimates  (22)  and  (25)  can  be  used  to  evaluate 
the  errors  | qt  —  qi,h\  if  the  remainder  terms  IZi  and  7 Zi  are  “sufficiently”  small; 
eventually,  they  can  be  used  to  correct  the  output  functionals  Si^h  as  =  qi,h+£qi 
or  qith  =  q%,h  +  £qi  ■ 

4  Error  estimators  with  B— spline  basis 

In  this  Section,  following  from  the  goal-oriented  a  posteriori  error  estimates  recalled  in 
Sec. 3,  we  propose  the  error  estimators  for  the  Galerkin  approximation  method  with  the 
use  of  a  B-spline  basis. 

4.1  Numerical  approximation:  B— spline  basis 

We  use  univariate  B-splinc  basis  for  the  construction  of  the  approximation  space  14,; 
we  refer  the  reader  to  [38]  and  also  to  [15]  for  a  general  overview  of  the  definition  and 
construction  of  B-spline  basis  as  well  as  their  properties.  We  consider  B-spline  basis 
defined  in  a  parametric  domain,  hereafter  denoted  by  f 2  =  (0, 1),  with  associated  knot 
vectors  S  =  {£,j}JLl7  with  the  jth  knot,  for  some  m  =  n  +  p  +  1,  where  n  is  the 
total  number  of  degrees  of  freedom  and  p  is  the  polynomial  order  under  consideration. 
Specifically,  we  consider  open  knots  vectors  S  with  internal  equally  spaced  Ne  knot 
spans,  for  which  the  total  number  of  degrees  of  freedom  is  n  =  Ne  +  p;  Ne  can  be 
interpreted  as  the  number  of  “elements”  e  =  1, . . . ,  Ne  of  size  h  in  which  12  is  partitioned. 
In  this  case,  5  =  Eh.p  :=  {{0}p+1, . . . ,  . . . ,  {1}P+1}  with  =  (j  -  p  -  1  )h  for 

j  =  p  +  2, . . . ,  n  with  h  :=  1  /Ne. 

The  jth  B-spline  basis  function,  BhtPlj(£),  is  defined  in  12  by  using  the  knot  vector  S 
with  the  Cox-de  Boor  recursion  formula  (see  [38]  and  [15]);  the  subscripts  h  and  p  in 
the  basis  functions  refer  to  the  number  of  internal  knot  spans  (Ne  “elements” )  and  their 
polynomial  order  p.  Examples  of  such  basis  functions  for  p  =  1,  2,  3  are  highlighted  in 
Fig.2(top-right),  (bottom-left)  and  (bottom-right),  respectively,  in  the  case  Ne  =  5. 
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Figure  2:  B-spline  basis  {Bh,P,i{f,)}Ij=\P  (left)  and  enhanced  basis  {B2h,p,j(Q}^=i+P 
(center)  and  {Bh,P+ij(£,)}i^iP+1  (right)  for  p  =  1  (top)  andp  =  2  (bottom)  with  £  £  f2; 
Ne  =  5  is  considered. 


We  observe  that  globally  high-order  continuous  B-spline  basis  functions  with  compact 
support  can  be  easily  constructed  by  introducing  a  relatively  small  number  of  additional 
degrees  of  freedom;  indeed,  n  =  Ne  +  p  with  respect  to  the  classic  Finite  Element 
method  for  which  C°  globally  continuous  basis  yields  n  =  pNe  +  1  degrees  of  freedom. 
Specifically,  the  basis  functions  so  far  constructed  are  Cp_1  globally  continuous  and  C°° 
in  each  “element” .  Additionally,  we  have  that  the  support  of  each  basis  function  is  p  + 1 
knot  spans  (“elements”). 

Remark  4.1  For  a  given  knot  vector 2,  the  B-spline  basis  of  order  p  =  1  ({Bh.i.jif)}^!^1 ) 
coincides  with  the  classic  linear  Finite  Element  Lagrangian  basis. 

The  geometrical  map  from  the  parametric  domain  f l  to  the  physical  domain  f 2, 
s(£)  :  — >  ft,  is  defined  by  using  the  B-spline  basis  {Bh,p,j(0}]j= l  and  n  control 

n 

points  Q  €  R  as  s(£)  :=  ^  Cj-  We  assume  that  Q+i  >  (j  for  j  =  1  ,...,n 

3= 1 

in  order  to  have  an  invertible  map,  with  (fi  =  0  and  =  L\  it  is  convenient  for  this 
one-dimensional  problem  to  distribute  the  control  points  Q  such  that  a  linear  map  is 
obtained. 

The  functional  space  Xh,P  '■=  span({B^iPJ  (£)}j=i )  °f  dimension  rih,p  =  Ne  +  p  is 
originated  by  the  B-spline  basis  for  £  £  f l,  for  which  a  generic  function  Vh,P(f)  £  XhtP 
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reads  Vh,P(£)  =  ^^JBh,P,j{C)vj  f°r  some  control  variables  Vj  £  R.  The  function  Vh,P(£) 
i= i 

mapped  in  the  physical  domain  Q  reads  Vh,P{s)  =  7a,p(£)os(0  •  Due  to  the  invertibility 

of  the  geometrical  map,  the  functions  Vh,p  and  Vh,P  will  be  indifferently  interchanged  as 
well  for  the  space  X h,P  with  its  correspondent  Xh,p  in  Q. 

It  follows  that  the  Galerkin  approximated  solutions  of  the  primal  and  dual  problems 
(Eqs.(ll)  and  (18)),  which  we  will  indicate  as  Oh, pi  zi,h,P  in  order  to  highlight  their 
dependency  on  the  “element”  size  h  and  order  p,  belongs  to  the  functional  space 
Vh.p  '■=  £  Xh,p  :  Vh,P |r  =  o|  C  V  of  dimension  nvh  p  <  rih,P ;  note  that  with  this 

new  notation,  the  space  Vh  =  Vh,P  and  the  output  functionals  =  qi,h,P- 
For  a  more  comprehensive  analysis  of  the  use  of  B-splines  basis  (eventually  NURBS 
basis)  and  the  related  advantages  in  the  context  of  analysis,  we  refer  the  reader  to  [7] 
and  [15]  (and  to  the  references  therein  indicated). 

We  solve  the  nonlinear  primal  problem  (14)  by  means  of  the  Newton-Raphson  method 
[39].  In  particular,  the  approximated  solution  0h,P  £  Vh,P  is  achieved  by  means  of 
a  converging  sequence  of  intermediate  solutions  0k  £  Vh,P,  obtained  by  solving  for 
k  =  0, 1, . . .  the  following  tangent  (linear)  problem  with  the  consequent  update: 

find  S0lp  £  Vh,p  :  a\elp){50lp,  <j>h>p)  =  -Rpr(0kh,p)(cfihtP)  V<f>h,p  £  Vh,p, 

(28) 

okh^=^,P  +  K,P, 

where  0°>p  £  VKp  is  a  prescribed  initial  guess  of  the  solution  “sufficiently”  close  to  the 
solution  0h,p'i  the  bilinear  form  a'{0){ •,  •)  is  given  in  Eq.(18),  while  the  primal  residual 
Rpr(-)(-)  in  Eq.(12).  The  satisfaction  of  the  stopping  criterium  ||Rfe+1||  <  fo(Ar_R||R°|| 
terminates  the  procedure  for  some  tol^R  -C  1,  where  Rfe  :=  {Rpr  (0k,P)(Bh,P,j)}j=i’P  ■ 

4.2  Error  estimators 

We  propose  and  consider  the  error  estimators  for  the  output  functionals  qi  for  i  = 
1, . . . ,  IV  by  using  the  a  posteriori  error  estimates  (22)  and  (25)  with  the  approximation 
method  presented  in  the  previous  Section.  Due  to  the  difficulties  in  properly  evaluating 
the  remainder  terms  IZi  and  IZi  of  Eqs.(24)  and  (27),  we  derive  error  estimators  by  using 
only  the  terms  Si  and  Si  of  Eqs.(23)  and  (26). 

We  remark  that  while  the  terms  Si  =  Si(6,Zi )  depend  on  both  the  exact  primal  and 
dual  solutions,  Si  =  Si(zi)  depends  only  on  the  exact  dual  solutions.  In  both  cases, 
in  order  to  derive  and  evaluate  the  error  estimators,  we  need  computable  variables  in 
place  of  the  exact  primal  and  dual  solutions.  The  most  natural  choice  for  replacing 
the  exact  solutions  with  the  approximated  ones  0h,P  and  Zi,h,P  £  Vh,P  is,  however,  not 
effective,  since  Si(0h,  zi,h,P)  =  Si(zi,h,P)  =  0  due  to  the  Galerkin  orthogonality  property 
of  the  primal  and  dual  problems  (11)  and  (18)  (see  also  [5]).  For  this  reason  the  exact 
solutions  need  to  be  replaced  by  computable  enhanced  approximated  solutions;  see  e.g. 
[8]  for  a  wider  discussion1 . 

A  first  couple  of  error  estimators  is  obtained  by  considering  enhanced  dual  solutions 
for  the  terms  Si(zi).  A  first,  standard  choice,  consists  in  building  the  enhanced  functional 
space  V2 h,P  by  uniformly  refining  the  knot  spans  internal  to  S  such  that  the  dimension 

An  alternative  possibility  consists  in  using  approximated  solutions  (or  gradients  of  the  solutions) 
recovered  from  Of,  p  and  Zij,p  S  V/f  /, ;  see  e.g.  [10]. 
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of  the  enhanced  space  X2h,p  is  U2h,P  =  2 Ne+p  =  2 nh,P  —  p;  this  procedure  is  called  knot 
insertion  ([38])  and  we  can  regard  it  as  “mesh”  refinement.  However,  knot  insertion 
does  not  coincide  with  the  h-refinement  of  the  Finite  Element  method  of  order  p ,  since 
the  new  basis  obtained  is  Cp~l  continuous  through  the  newly  inserted  knots  and  not 
C°;  they  do  coincide  only  for  the  case  p=  1. 

As  alternative,  we  propose  to  use  a  higher  order  enhanced  dual  solution  with  respect  to 
the  order  p  considered  for  the  primal  solution  9htP  G  V/j,p.  With  this  aim,  in  order  to 
build  the  higher  order  enhanced  space  Vh,p+i,  we  consider  the  order  elevation  of  the  orig¬ 
inal  B-spline  basis,  but  without  increasing  the  multiplicity  of  the  internal  knots  of  the 
original  knot  vector  ([38]);  only  the  multiplicity  of  the  external  knots  {0}  and  {1}  is 
increased  by  one  to  obtain  the  new  knot  vector  S/jiP+ 1  :=  {{0}p+2, . . .  . . . ,  {1}P+2}. 

The  result  is  a  B-spline  basis  with  increased  Cp  continuity  through  the  internal  knots 
£j  of  the  knot  vector  EhtP+i2-  This  procedure  corresponds  to  a  step  of  the  “pure”  k- 
refinement  approach  introduced  for  Isogeometric  Analysis  [15];  however,  this  can  not  be 
identified  with  the  standard  p-refinement  of  the  Finite  Element  method  which  would 
maintain  C°  continuity  through  the  knots.  Finally,  the  higher  order  enhanced  space 

reads  Vh.p+i  :=  ^Vh,P+i  G  XhtP+i  ■  f/j,P+i|r  =  oj,  where  the  space  Xh,P+ i  assumes 
dimension  rih,p+ i  =  Ne  +  p  +  1  =  +  1. 

Examples  of  the  so  obtained  enhanced  B-spline  basis  are  displayed  in  Fig. 2  for  the  cases 
p  =  1,2  and  Ne  =  5. 

The  enhanced  dual  solutions  Zi,2h,P  G  V2 h,P  and  Zi,h,P+ 1  G  Vh,P+i  are  obtained  by  solving 
the  following  dual  problems,  respectively: 

find  Zi^2h,P  G  V2/i,P  :  a  ($/i,p)(£/,2/i,pi  02/i,p)  ~  ^(^/i,p)(^2/i,p)  V02/i,p  G  V2 /i,P,  (29) 

find  Zi^h.p- t-i  G  Vh,P+\  :  a  ($/i,p)(z/,/i,p+i;  </>/&, p-t-i)  —  li{9h,P){.4>h,P+ 1)  ^tyh^p+i  G  V/i,p-j-i, 

(30) 

for  i  =  1, . . . ,  N,  where  9h,P  G  Vh,P  is  the  approximated  solution  of  the  primal  problem 
(not  enhanced). 

On  this  basis,  from  Eq.(26),  we  define  the  error  estimators  A i^h,P  and  as, 

respectively: 

A i,2h,P '■=  £i{zi,2h,P)  i  =  l,...,N,  (31) 

^i,h,p+l  :=:  'l  =  1,  •  ■  •  ,  Al,  (32) 

where  the  ith  index  refer  to  the  output  functional  qt.  The  error  estimator  A it2h,P  co¬ 
incides  with  the  one  considered  in  [46],  proposed  for  two-dimensional  B-splines  and 
T-splines. 

We  provide  two  additional  error  estimators  from  the  term  £i(9,Zi)  of  Eq.(23),  for 
which  the  concept  of  enhanced  dual  solutions  introduced  so  far  can  be  extended  to 
the  primal  solution  G  Vh,P-  However,  we  observe  that  in  order  to  obtain  such 
enhanced  solutions  6*2/1, P  G  V2 h,P  and  9hlP+ 1  G  Vh,P+i,  nonlinear  problems  need  to  be 
solved  by  means  of  the  Newton-Raphson  method.  In  order  to  avoid  this  inefficiency  in 
the  evaluation  of  the  errors  on  the  outputs  qi  associated  with  the  current  primal  solution 
9h,P  G  Vh,P ,  we  propose  an  alternative  approach  based  on  the  tangent  problem  (28)  of 
the  Newton-Raphson  method.  In  particular,  we  obtain  the  enhanced  primal  solutions 
6*2/1, P  G  V2 h,P  and  6*/j,p+i  G  Vh,p+\  by  solving  the  following  tangent  (linear)  problems, 

2  In  the  standard  order  elevation  procedure  the  multiplicity  of  the  internal  knots  is  increased  in  order 
to  preserve  the  discontinuity  of  the  derivatives  of  the  enhanced  basis  with  respect  to  the  original  one. 
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respectively: 


filld  502h,p  ^  ^2 h,p  •  G  {Oh,p)  (S@2h,p  5  4)2h,p)  —  [Ph.p'){t4>2h,p]  ^4*2 h,p  £  V 2h.p , 


@2 h,p  ~  II 2^p0h,p  +  S02h,p, 


(33) 


find  (z  Vfr)P+l  .  G  (@h,p)(.3@h,p+li  4*h,p-\-l)  —  R^  (@h,p)((fih,p+ 1)  '^4)h,p+ 1  f=  V/^p-j-1, 


0h,P+1  =  n^+10fciP  +  <j0h,p+i, 

(34) 

for  the  given  solution  of  the  primal  problem  6htP  €  Vh,p]  nh  ^  :  V/ijP  — >•  V2h,,P  and 

II h’p+1  ■  Vh,p  — t  V/i,p+i  represent  projection  operators  of  functions  Vh,P  S  Vh.p  into  the 
spaces  V2/t,P  and  V/liP+i,  respectively. 

By  recalling  Eq.(23)  and  the  enhanced  dual  solutions  of  Eqs.(29)  and  (30),  we  define 
the  error  estimators  Aj;2/i,P  and  Aith:P+i  as: 

A^2/l,P  ■ —  ^(^2/l,P;  ^Z,2/l,P)  ^  —  !)•••>  A)  (35) 

Ai^h.p- j-1  :—  &i{@h,p+ 1;  ^2,/i,P-}-l)  i  ~  1,  .  .  .  ,  iV .  (36) 


We  remark  that  the  error  estimators  can  be  conveniently  and  immediately  adopted 
in  the  case  of  linear  problems. 


Remark  4.2  Following  from  Remark  J^.l,  the  error  estimators  A^2/i,p  and  A ^2/1, P  for 
p  =  1  can  be  directly  used  to  estimate  errors  associated  with  the  classic  Finite  Element 
method  with  linear  basis. 


We  note  that  the  error  estimators  do  not  represent  rigorous  error  bounds  for  the 
errors  on  the  output  functionals  \qi  —  qi,h\,  since  the  remainder  terms  IZi  and  IZi  are  not 
included  in  the  estimators  and  the  exact  primal  and  dual  solutions  are  replaced  with 
enhanced  solutions.  In  order  to  evaluate  the  capability  of  the  error  estimators  (31), 
(32),  (35)  and  (36)  to  bound  and  provide  indications  of  the  error  (21),  we  introduce  the 
following  effectivity  indexes: 


Pi,2h,p 

Vi,2h,p 


Ai,2h,p 
| qi  qi.h.p  | 

A  i,2h,p 

|  qi  qi,h,p  | 


Vi,h,p+ 1 
Vi,h,p+ 1 


Ai^iP-(-l 
|  qi  Qi,h,p  | 
Aj^p+i 
|  qi  Qi,h,p  | 


(37) 


for  i  =  1, ...  ,N.  Effectivity  indexes  larger  than  1  indicate  the  capability  of  the  error  es¬ 
timator  to  bound  the  error;  however,  sharp  estimators  should  exhibit  effectivity  indexes 
close  to  1.  Asymptotic  effectivity  indexes  are  obtained  for  large  number  of  degrees  of 
freedom  n\>h  . 


Remark  4.3  Although  we  are  only  dealing  with  a  one-dimensional  problem ,  it  is  still 
interesting  to  discuss  the  computational  costs  associated  with  the  evaluation  of  the  es¬ 
timators.  With  this  respect,  the  estimator  Aj^p+i  is  the  most  convenient  one,  since 
the  primal  residual  is  evaluated  with  the  enhanced  dual  solutions  Zi^h,p+i  for  which  only 
an  additional  degree  of  freedom  is  added  with  respect  to  0h,p-  Conversely,  the  estimator 
A i,2h,p  is  the  most  expensive  to  evaluate  since  both  the  primal  and  dual  residuals  need 
to  be  evaluated  with  the  enhanced  primal  and  solutions  02h,p  ond  Zit2h}P,  each  containing 


11 


Test 

El 

L 

px 

Py 

M 

i.i 

10.0 

1.00 

0.0 

100 

0.0 

1.2 

” 

” 

1  n2EI 

2  \T? 

10.0 

” 

1.3 

11 

” 

3  7 ?EI 

2  41? 

” 

” 

2 

3.60 -10^ 

7T 

-2.062648 

0.0 

0.0 

3.1 

0.50 

4.00 

0.150 

-1.00 

1  7T  El 

~To  j, 

3.2 

0.50 

4.00 

1.00 

0.0 

0.0 

4 

10.0 

(27r  —  a)R 

50.0 

0.0 

0.0 

Table  1:  Data  for  the  Test  problems  1.1  1.3,  2,  3. 1-3. 2  and  4. 


Test 

So  (s) 

“Exact” 

Ne  =  2 

4 

8 

16 

32 

64 

128 

1.1 

7 r  s 

5 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

1.2 

ii 

5 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

1.3 

ii 

4 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

2 

—7 r  s 

4 

4,  4 

4,  4 

4,  4 

4,  4 

4,  4 

4,  4 

4,  4 

3.1 

— 7T  S 

7 

10,  16 

11,  13 

10,  6 

7,  7 

7,  7 

7,  7 

7,  7 

3.2 

■k/2  s 

5 

4,  4 

4,  4 

5,  5 

5,  5 

5,  5 

5,  5 

5,  5 

4 

—  cos(tts) 

5 

4,  6 

5,  7 

5,  6 

5,  5 

5,  5 

5,  5 

5,  5 

Table  2:  Number  of  Newton-Raphson  iterative  steps  starting  from  the  initial  solution 
90(s)  for  the  “exact11  solution  and  the  approximate  ones  with  Ne  £  {2r}£=2  ( p  =  1, 
p  =  2)  for  the  Test  problems  1. 1-1.3,  2,  3.1  -3.2  and  4;  the  pairs  in  the  entries  of  the 
table  indicate  the  steps  for  p  =  1  and  p  =  2. 


a  number  of  degrees  of  freedom  about  (2 Ne  +p)/(Ne  +  p)  times  higher.  Among  these, 
intermediate  costs  are  obtained  with  Ai^h,P  and  Aith,P+i-  Indeed,  A it2h,P  only  require 
the  evaluation  of  the  primal  residual  but  with  the  enhanced  dual  solution  Zi^h.p,  while 
Ai:h,p+i  evaluates  both  the  primal  and  dual  residuals  with  the  enhanced  primal  and  dual 
solutions  Oh.p+i  and  Zith,P+ 1-  In  general,  the  cost  associated  with  Aj,^h,P  could  be  higher 
than  the  one  for  Ai^,p+i  for  a  sufficiently  large  Ne. 

Additionally,  the  cost  associated  with  the  computations  of  the  enhanced  primal  and  dual 
solutions  should  be  taken  into  account.  With  this  respect  the  estimators  Ai2h,P  and 
Ai,h,p+i  only  require  the  solution  of  the  linear  dual  problems  (29)  and  (30),  while  for 
A it2h,P  and  Aith,p+i,  also  the  enhanced  primal  ones  (33)  and  (34)  need  to  be  solved, 
with  consequent  increased  computational  costs. 

5  Results  and  discussion 

In  this  Section  we  report  the  numerical  results  obtained  using  the  proposed  error  es¬ 
timators  applied  to  the  solution  of  various  Test  problems  of  the  same  type  as  the  one 
introduced  in  Sec. 2.  The  error  estimators  are  evaluated  and  analyzed  in  terms  of  their 
effectivity  indexes;  a  discussion  follows  from  these  results. 
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0 


0 


Figure  3:  Undeformed  configuration  (dashed  line, - )  and  deformed  ones  (continu¬ 

ous  lines,  — )  for  Tests  1.1- 1.3  (top-left),  2  (top-right),  3. 1-3. 2  (bottom-left)  and  4 
(bottom-right). 


5.1  Numerical  results 

We  study  four  Test  problems,  Tests  1-4;  Tests  1.1  1.3  and  3. 1-3. 2  represent  the  same 
problems  but  with  different  load  conditions.  The  data  used  for  the  various  Test  problems 
are  reported  in  Table  1  (see  Sec. 2  for  the  definition  of  the  notation).  We  note  that  Test  4 
does  not  correspond  to  a  cantilever  problem,  but  to  a  simply  supported  rod  problem. 
However,  as  both  Py  and  M  are  set  to  zero,  we  can  make  use  of  the  equations  for  the 
cantilever  problem  with  no  changes;  in  particular,  for  this  Test  we  choose  a  =  4-7r/5  and 
R  =  1.  The  “exact”  solutions  of  the  Test  problems  are  plotted  in  Fig. 3  in  terms  of  the 
displacements  u(s)  and  v(s).  Such  solutions  are  obtained  using  the  Newton- Raphson 
method  with  initial  solutions  as  the  ones  indicated  with  9o(s)  in  Table  2. 

We  solve  the  Test  problems  using  B-spline  bases  of  orders  p  =  1  and  p  =  2  for 
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Figure  4: 
derivatives 
lines, - ) 


Test  3.2.  Comparison  of  the  approximated  solutions  0h,P  (top)  and  their 
dOh 

— (bottom)  with  the  “exact”  ones  (continuous  lines,  — )  forp  =  1  (dashed 
as 

and  p  =  2  (dash-dotted  lines,  —  •);  Ne  =  4  (left)  and  Ne  =  8  (right). 


decreasing  values  of  h  =  1  /Ne-  in  particular,  we  assume  Ne  =  2, 4,  8, . . . ,  128  (Ne  £ 
{2r}J=2).  A  7-point  Gauss  quadrature  rule  is  used  for  the  evaluation  of  the  integrals 
in  each  element  e  =  1  ,...,_ZVe.  Also,  the  tolerance  for  the  stopping  criterion  of  the 
Newton-Raphson  method  (28)  is  set  to  tolNR  =  10-12. 

Additionally,  we  assume  the  approximated  solutions  obtained  by  solving  the  primal 
problem  with  a  B-spline  basis  of  order  p  =  4  with  an  open  knot  vector  S  partitioned 
into  Ne  =  1,024  equally  spaced  “elements”  as  the  “exact”  solutions. 

In  Fig. 4  we  compare  the  approximated  solutions  and  their  derivatives,  obtained  for 
Test  3.2  with  p  =  1,2  and  Ne  =  4,  8,  against  the  “exact”  solutions. 

In  Table  2  we  report  the  number  of  Newton-Raphson  iterative  steps  required  for  conver¬ 
gence  to  the  approximate  primal  solutions  9h,P  £  Vh,P  hr  comparison  to  those  obtained 
for  convergence  to  the  “exact”  solution.  We  can  observe  that  the  obtained  number  of 
steps  is  equal  to  the  one  obtained  for  convergence  to  the  “exact”  solution  in  most  of 
the  cases,  with  the  exception  of  Tests  3  and  4  when  a  small  number  of  “elements”  Ne 
is  considered. 

In  Tables  3-9  we  report  the  effectivity  indexes  (37)  of  the  error  estimators  associated 
with  all  the  Test  problems  and  output  functionals  qi  for  i  =  1, . . . ,  N  =  4  (the  index 
i  is  dropped  in  the  effectivity  indexes,  since  deducible  from  the  context)  for  all  the 
approximations  Vh,P  obtained  with  Ne  £  {2r}J=2  and  p  =  1,2.  The  values  of  the 
“exact”  output  functionals  qi  are  reported  as  well  as  the  asymptotic  convergence  orders 
of  the  errors  \qi  —  qi,h,P\  with  respect  to  h ;  these  correspond  to  about  2  and  4  in  the 
cases  p  =  1  and  p  =  2,  respectively.  The  output  functional  q 4  is  not  discussed  for  Test  4 
since  in  this  case  54  =  0  and  qith,P  =  0  except  for  round-off  errors. 
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P  =  1 


p  =  2 


Ne 

Iqi  Ql,h,p\ 

V2h,p 

Vh,p+ 1 

92  h,p 

T]h,p+1 

2 

1.3257e  +  00 

0.79 

1.18 

0.79 

0.60 

4 

3.2983e  —  01 

0.89 

2.05 

0.89 

1.38 

8 

8.0124e  —  02 

0.95 

2.68 

0.95 

1.83 

16 

1.9863e  —  02 

0.97 

3.00 

0.97 

2.03 

32 

4.9551e  —  03 

0.97 

3.15 

0.97 

2.12 

64 

1.2381e  —  03 

0.97 

3.22 

0.97 

2.16 

128 

3.0949e  -  04 

0.98 

3.26 

0.98 

2.18 

Conv.  Order  =  2.00 


Ne 

\l2  -  l2,h,p\ 

V2h,p 

Wh,p+1 

V2h,p 

Vh,p+ 1 

2 

4.5317e  —  02 

0.75 

0.82 

0.75 

0.93 

4 

1.0614e  —  02 

0.75 

0.96 

0.75 

1.11 

8 

2.6207e  —  03 

0.75 

0.99 

0.75 

1.15 

16 

6.5324e  —  04 

0.75 

1.00 

0.75 

1.16 

32 

1.6319e  —  04 

0.75 

1.00 

0.75 

1.17 

64 

4.0791e  —  05 

0.75 

1.00 

0.75 

1.17 

128 

1.0197e  —  05 

0.75 

1.00 

0.75 

1.17 

Conv.  Order  =  2.00 


Ne 

\l3  ~  93, ft, pi 

V2h,p 

Vh,p+ 1 

V2h,p 

Vh,p+1 

2 

2.7692e  —  02 

0.72 

0.98 

0.72 

0.68 

4 

6.3033e  —  03 

0.74 

1.05 

0.74 

0.50 

8 

1.5299e  —  03 

0.75 

1.02 

0.75 

0.41 

16 

3.7993e  —  04 

0.75 

1.01 

0.75 

0.37 

32 

9.4829e  —  05 

0.75 

1.00 

0.75 

0.36 

64 

2.3698e  —  05 

0.75 

1.00 

0.75 

0.36 

128 

5.9239e  —  06 

0.75 

1.00 

0.75 

0.36 

Conv.  Order  =  2.00 


Ne 

\Q4  Q4,h,p\ 

V2h,p 

r)h,p+ 1 

V2h,p 

Vh,p+1 

2 

4.3420e  —  02 

0.81 

0.88 

0.81 

0.75 

4 

1.0687e  —  02 

0.77 

0.97 

0.77 

0.67 

8 

2.6301e  —  03 

0.75 

1.00 

0.75 

0.64 

16 

6.5463e  —  04 

0.75 

1.00 

0.75 

0.62 

32 

1.6348e  —  04 

0.75 

1.00 

0.75 

0.62 

64 

4.0858e  -  05 

0.75 

1.00 

0.75 

0.61 

128 

1.0214e  —  05 

0.75 

1.00 

0.75 

0.61 

Conv.  Order  =  2.00 


9i  =  17.954,  92  =  1.4303, 


Ne 

|#1  Ql,h,p\ 

V2  h,p 

f]h,p+l 

92/1,1) 

Vh,p+ 1 

2 

1.3831e  —  01 

0.51 

0.97 

0.51 

1.73 

4 

5.4334e  -  03 

0.42 

1.65 

0.42 

6.83 

8 

3.1742e  -  04 

0.40 

1.63 

0.40 

8.32 

16 

1.9405e  -  05 

0.39 

1.60 

0.39 

8.04 

32 

1.1957e—  06 

0.39 

1.59 

0.39 

7.70 

64 

7.4366e  -  08 

0.39 

1.58 

0.39 

7.49 

128 

4.6413e  —  09 

0.39 

1.58 

0.39 

7.38 

Conv.  Order  =  4.00 


Ne 

\l2  -  ?2,h,p| 

T]2h,p 

IjhyP+l 

V2h,p 

Vh,p+ 1 

2 

9.0134e  -  04 

0.92 

0.66 

0.92 

0.18 

4 

7.9694e  -  05 

0.93 

0.96 

0.93 

0.83 

8 

5.3464e  -  06 

0.94 

0.99 

0.94 

1.01 

16 

3.3891e  -  07 

0.94 

1.00 

0.94 

0.87 

32 

2.1228e  —  08 

0.94 

1.00 

0.94 

0.75 

64 

1.3272e  —  09 

0.94 

1.00 

0.94 

0.69 

128 

8.2955e  —  11 

0.94 

1.00 

0.94 

0.66 

Conv.  Order  =  4.00 


Ne 

1*73  -  Q3,h,p\ 

H2h,p 

nh,p+ 1 

H2h,p 

Vh,p+ 1 

2 

2.0920e  -  05 

1.64 

9.74 

1.64 

62.5 

4 

8.2291e  —  06 

1.02 

0.45 

1.02 

39.6 

8 

4.5144e  -  08 

0.92 

3.78 

0.92 

526 

16 

6.1625e  —  09 

0.91 

0.39 

0.91 

228 

32 

5.7855e  —  10 

0.93 

0.89 

0.93 

143 

64 

4.0226e  —  11 

0.94 

0.97 

0.94 

124 

128 

2.5867e  -  12 

0.94 

0.99 

0.93 

119 

Conv.  Order  =  3.96 


Ne 

\Q4  Q4,h,p  | 

n2h,p 

Vh,p+ 1 

H2h,p 

Vh,p+ 1 

2 

2.0822e  -  03 

0.98 

0.93 

0.98 

0.28 

4 

8.4967e  -  05 

0.95 

1.08 

0.95 

2.82 

8 

5.0163e  -  06 

0.94 

1.03 

0.94 

3.87 

16 

3.0663e  -  07 

0.94 

1.01 

0.94 

3.72 

32 

1.8916e  —  08 

0.94 

1.00 

0.94 

3.51 

64 

1.1771e  —  09 

0.94 

1.00 

0.94 

3.37 

128 

7.3478e  -  11 

0.94 

1.00 

0.94 

3.29 

Conv.  Order  =  4.00 


93  =  -0.55500,  94  =  0.81061 


Table  3:  Test  1.1.  Errors  and  effectivity  indexes  772/1, p,  Vh,p+i,  172/1, p  and  rjh,p+i  for  the 
output  functionals  gi,  g2,  g3  and  g4  with  Ne  £  {2r}J_2  and  p  =  1,2;  the  asymptotic 
convergence  orders  and  “exact”  values  of  the  output  functionals  are  reported. 
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P  =  1 


p  =  2 


Ne 

Iqi  Ql,h,p\ 

V2h,p 

Vh,p+1 

ri2h,p 

T]h,p+1 

2 

2.9640e  —  01 

0.22 

1.39 

0.22 

0.73 

4 

7.3836e  -  02 

0.20 

1.55 

0.20 

0.83 

8 

1.8444e  —  02 

0.20 

1.64 

0.20 

0.88 

16 

4.6100e  —  03 

0.20 

1.69 

0.20 

0.91 

32 

1.1524e  —  03 

0.20 

1.71 

0.20 

0.92 

64 

2.8810e  —  04 

0.20 

1.72 

0.20 

0.93 

128 

7.2026e  —  05 

0.20 

1.73 

0.20 

0.93 

Conv.  Order  =  2.00 


Ne 

|92  -  92,h,pl 

V2h,p 

Wh,p+1 

V2h,p 

Vh,p+1 

2 

8.0199e  —  03 

0.82 

1.15 

0.82 

0.68 

4 

1.8738e  —  03 

0.77 

1.04 

0.77 

0.48 

8 

4.6054e  —  04 

0.75 

1.01 

0.75 

0.34 

16 

1.1464e  —  04 

0.75 

1.00 

0.75 

0.27 

32 

2.8631e  —  05 

0.75 

1.00 

0.75 

0.23 

64 

7.1557e  —  06 

0.75 

1.00 

0.75 

0.21 

128 

1.7888e  —  06 

0.75 

1.00 

0.75 

0.20 

Conv.  Order  =  2.00 


Ne 

l?3  -  Q3.h,p\ 

V2h.p 

r)h,p+ 1 

V2h,p 

Vh,p+1 

2 

1.7068e  —  02 

0.73 

1.03 

0.73 

0.65 

4 

4.3594e  -  03 

0.75 

1.01 

0.75 

0.71 

8 

1.0957e  —  03 

0.75 

1.00 

0.75 

0.74 

16 

2.7430e  —  04 

0.75 

1.00 

0.75 

0.76 

32 

6.8597e  —  05 

0.75 

1.00 

0.75 

0.77 

64 

1.7151e  —  05 

0.75 

1.00 

0.75 

0.77 

128 

4.2878e  —  06 

0.75 

1.00 

0.75 

0.77 

Conv.  Order  =  2.00 


Ne 

\Q4  Q4,h,p\ 

V2h,p 

r)h,p+ 1 

V2 h.p 

Vh,p+1 

2 

3.0014e  —  02 

0.77 

1.03 

0.77 

0.64 

4 

7.4033e  -  03 

0.75 

1.01 

0.75 

0.70 

8 

1.8446e  —  03 

0.75 

1.00 

0.75 

0.73 

16 

4.6075e  -  04 

0.75 

1.00 

0.75 

0.75 

32 

1.1516e  —  04 

0.75 

1.00 

0.75 

0.75 

64 

2.8789e  —  05 

0.75 

1.00 

0.75 

0.76 

128 

7.1973e  —  06 

0.75 

1.00 

0.75 

0.76 

Conv.  Order  =  2.00 


gi  =  3.7685,  ®  =  0.75854, 


Ne 

|#1  Ql,h,p\ 

V2h,p 

f]h,p+l 

V2h,p 

Vh,p+ 1 

2 

6.7273e  -  04 

0.54 

0.72 

0.54 

1.82 

4 

4.5086e  -  05 

0.58 

1.21 

0.58 

7.38 

8 

2.6654e  -  06 

0.57 

1.33 

0.57 

16.7 

16 

1.6335e  —  07 

0.57 

1.37 

0.57 

34.5 

32 

1.0153e  —  08 

0.57 

1.38 

0.57 

69.9 

64 

6.3364e  -  10 

0.57 

1.39 

0.57 

140 

128 

Conv.  Order  =  4.12 


Ne 

l?2  ~  92, h, p] 

V2h,p 

IjhyP+l 

V2h,p 

Vh,p+1 

2 

3.6097e  -  05 

0.94 

0.76 

0.94 

11.4 

4 

2.3489e  -  06 

0.94 

0.98 

0.94 

60.2 

8 

1.3765e  -  07 

0.94 

1.00 

0.94 

146 

16 

8.4227e  -  09 

0.94 

1.00 

0.94 

311 

32 

5.2336e  —  10 

0.94 

1.00 

0.94 

639 

64 

3.2659e  —  11 

0.94 

1.00 

0.94 

1293 

128 

Conv.  Order  =  4.00 


Ne 

193  -  93, A, pi 

H2h,p 

nh,p+ 1 

H2h,p 

Vh,p+ 1 

2 

7.4043e  -  05 

0.94 

0.92 

0.94 

3.26 

4 

4.4798e  -  06 

0.94 

1.00 

0.94 

15.1 

8 

2.6965e  —  07 

0.94 

1.00 

0.94 

34.7 

16 

1.6661e  —  08 

0.94 

1.00 

0.94 

72.0 

32 

1.0381e  —  09 

0.94 

1.00 

0.94 

147 

64 

6.4829e  —  11 

0.94 

1.00 

0.94 

297 

128 

Conv.  Order  =  4.00 


Ne 

\Q4  Q4,h,p\ 

r]2h,p 

Vh,p+ 1 

H2h,p 

Vh,p+ 1 

2 

3.8781e  -  06 

0.86 

1.55 

0.86 

79.0 

4 

7.1194e  —  07 

0.95 

0.77 

0.95 

158 

8 

3.8296e  -  08 

0.94 

0.95 

0.94 

420 

16 

2.2311e  —  09 

0.94 

0.99 

0.94 

942 

32 

1.3647e  -  10 

0.94 

1.00 

0.94 

1967 

64 

8.4783e  -  12 

0.94 

1.00 

0.94 

3998 

128 

Conv.  Order  =  4.01 


93  =  -0.14608,  94  =  0.46902 


Table  4:  Test  1.2.  Errors  and  effectivity  indexes  772/1, p,  Vh,p+i,  V2h,p  and  r)h,p+i  for  the 
output  functionals  q±,  (72,  93  and  <74  with  Ne  £  {2r}J_2  and  p  =  1,2;  the  asymptotic 
convergence  orders  and  “exact”  values  of  the  output  functionals  are  reported.  Errors 
and  effectivity  indexes  for  Ne  =  128  and  p  =  2  are  not  reported  since  roundoff  errors 
disrupt  the  convergence  order. 
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P  =  1 


p  =  2 


Ne 

Iqi  Ql,h,p\ 

V2h,p 

Vh,p+ 1 

V2h,p 

f]h,p+l 

2 

9.6686e  —  01 

0.09 

1.81 

0.09 

0.88 

4 

2.2514e  —  01 

0.21 

2.09 

0.21 

0.97 

8 

5.5864e  —  02 

0.23 

2.21 

0.23 

1.00 

16 

1.3945e  —  02 

0.24 

2.27 

0.24 

1.01 

32 

3.4849e  -  03 

0.24 

2.29 

0.24 

1.02 

64 

8.7115e  —  04 

0.24 

2.31 

0.24 

1.02 

128 

2.1778e  —  04 

0.24 

2.31 

0.24 

1.02 

Conv.  Order  =  2.00 


Ne 

\Q2  ~  92,A,p| 

V2h,p 

Vh,p+l 

V2h,p 

Vh,p+ 1 

2 

1.4640e  —  02 

0.54 

0.53 

0.54 

2.26 

4 

4.3532e  —  03 

0.71 

0.91 

0.71 

1.83 

8 

1.1161e  —  03 

0.74 

0.98 

0.74 

1.58 

16 

2.8063e  —  04 

0.75 

0.99 

0.75 

1.45 

32 

7.0257e  —  05 

0.75 

1.00 

0.75 

1.38 

64 

1.7570e  —  05 

0.75 

1.00 

0.75 

1.34 

128 

4.3929e  —  06 

0.75 

1.00 

0.75 

1.32 

Conv.  Order  =  2.00 


Ne 

\l3  ~  </3,A,p| 

V2h,p 

Vh,p+ 1 

V2h  ,p 

Vh,p+1 

2 

5.2081e  —  02 

0.77 

1.11 

0.77 

0.70 

4 

1.2768e  —  02 

0.75 

1.03 

0.75 

0.74 

8 

3.1931e  —  03 

0.75 

1.01 

0.75 

0.81 

16 

7.9846e  -  04 

0.75 

1.00 

0.75 

0.85 

32 

1.9963e  —  04 

0.75 

1.00 

0.75 

0.88 

64 

4.9908e  —  05 

0.75 

1.00 

0.75 

0.89 

128 

1.2477e  —  05 

0.75 

1.00 

0.75 

0.90 

Conv.  Order  =  2.00 


Ne 

\Q4  Q4,h,p\ 

V2h,p 

r)h,p+ 1 

V2 h,p 

Vh,p+1 

2 

2.1520e  —  02 

0.86 

1.03 

0.86 

0.62 

4 

4.6625e  —  03 

0.78 

1.01 

0.78 

0.68 

8 

1.1309e  —  03 

0.76 

1.00 

0.76 

0.76 

16 

2.8067e  —  04 

0.75 

1.00 

0.75 

0.81 

32 

7.0041e  -  05 

0.75 

1.00 

0.75 

0.84 

64 

1.7502e  —  05 

0.75 

1.00 

0.75 

0.86 

128 

4.3751e  —  06 

0.75 

1.00 

0.75 

0.86 

Conv.  Order  =  2.00 


9i  =  20.877,  92  =  1.7732, 


Ne 

|<?1  Ql,h,p\ 

92  h,p 

f]h,p+l 

V2h,p 

Vh,p+ 1 

2 

1.4633e  -  02 

0.21 

0.01 

0.21 

0.05 

4 

1.8363e  -  03 

0.55 

1.19 

0.55 

0.53 

8 

9.3996e  -  05 

0.51 

1.39 

0.51 

3.06 

16 

5.5185e  —  06 

0.50 

1.44 

0.50 

7.78 

32 

3.4011e  —  07 

0.49 

1.45 

0.49 

17.0 

64 

128 

Conv.  Order  =  4.02 


Ne 

\l2  ~  92, A, p| 

T]2h,p 

r)h,p+ 1 

V2h,p 

Vh,p+ 1 

2 

3.7772e  -  04 

0.93 

0.64 

0.93 

2.29 

4 

3.3988e  -  05 

0.95 

1.01 

0.95 

19.5 

8 

1.6319e  —  06 

0.94 

1.03 

0.94 

56.1 

16 

9.3868e  -  08 

0.94 

1.01 

0.94 

123 

32 

5.7736e  -  09 

0.94 

1.00 

0.94 

252 

64 

128 

Conv.  Order  =  4.02 


Ne 

1*73  -  93, A, pi 

H2h,p 

nh,p+ 1 

f]2h,p 

Vh,p+ 1 

2 

8.2262e  -  04 

0.91 

0.74 

0.91 

1.38 

4 

7.2521e  —  05 

0.95 

1.01 

0.95 

9.03 

8 

3.8584e  -  06 

0.94 

1.02 

0.94 

22.6 

16 

2.3013e  —  07 

0.94 

1.01 

0.94 

47.1 

32 

1.4232e  —  08 

0.94 

1.00 

0.94 

95.0 

64 

128 

Conv.  Order  =  4.02 


Ne 

\Q4  Q4,h,p\ 

'H2h,p 

Vh,p+1 

H2h,p 

Vh,p+ 1 

2 

5.1043e  -  04 

0.97 

1.15 

0.97 

0.53 

4 

8.1165e  —  06 

0.92 

1.66 

0.92 

31.6 

8 

6.0445e  -  07 

0.93 

1.09 

0.93 

62.1 

16 

4.1133e  —  08 

0.94 

1.02 

0.94 

117 

32 

2.6268e  -  09 

0.94 

1.00 

0.94 

233 

64 

128 

Conv.  Order  =  3.97 


94  =  —0.68420,  93  =  0.80464 


Table  5:  Test  1.3.  Errors  and  effectivity  indexes  772/1, p,  Vh,p+u  172/1, p  and  r)h,P+i  for  the 
output  functionals  gi,  g2,  g3  and  g4  with  Ne  £  {2r}J_2  and  p  =  1,2;  the  asymptotic 
convergence  orders  and  “exact”  values  of  the  output  functionals  are  reported.  Errors 
and  effectivity  indexes  for  Ne  =  64, 128  and  p  =  2  are  not  reported  since  roundoff  errors 
disrupt  the  convergence  order. 
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P  =  1 


p  =  2 


Ne 

Iqi  Ql,h,p\ 

V2h,p 

rlh,p+ 1 

V2h,p 

T]h,p+1 

2 

2.2404e  +  00 

1.84 

1.26 

1.84 

1.20 

4 

4.7824e  -  01 

2.30 

0.54 

2.30 

1.58 

8 

1.1172e  —  01 

2.48 

0.30 

2.48 

1.44 

16 

2.7423e  —  02 

2.53 

0.84 

2.53 

1.24 

32 

6.8241e  —  03 

2.55 

1.12 

2.55 

1.13 

64 

1.7040e  —  03 

2.55 

1.26 

2.55 

1.06 

128 

4.2589e  —  04 

2.55 

1.33 

2.55 

1.03 

Conv.  Order  =  2.00 


Ne 

|92  -  92,h,pl 

V2h,p 

Wh,p+1 

V2h,p 

Vh,p+ 1 

2 

1.0033e  —  01 

0.72 

0.89 

0.72 

0.89 

4 

2.7267e  —  02 

0.74 

0.98 

0.74 

0.97 

8 

6.8939e  —  03 

0.75 

1.00 

0.75 

0.99 

16 

1.7263e  —  03 

0.75 

1.00 

0.75 

1.00 

32 

4.3174e  —  04 

0.75 

1.00 

0.75 

1.00 

64 

1.0794e  —  04 

0.75 

1.00 

0.75 

1.00 

128 

2.6986e  —  05 

0.75 

1.00 

0.75 

1.00 

Conv.  Order  =  2.00 


Ne 

l?3  -  Q3,h,p\ 

V2h  ,p 

Vh,p+ 1 

V2h,p 

Vh,p+1 

2 

6.6589e  —  01 

0.79 

0.50 

0.79 

0.66 

4 

2.5881e  —  01 

0.76 

0.80 

0.76 

0.88 

8 

7.1712e  —  02 

0.75 

0.94 

0.75 

0.97 

16 

1.8363e  —  02 

0.75 

0.98 

0.75 

0.99 

32 

4.6178e  —  03 

0.75 

1.00 

0.75 

1.00 

64 

1.1561e  —  03 

0.75 

1.00 

0.75 

1.00 

128 

2.8913e  —  04 

0.75 

1.00 

0.75 

1.00 

Conv.  Order  =  2.00 


Ne 

1^4  Q&,h,p\ 

V2h,p 

Vh,p+ 1 

V2h,p 

Wh,p+1 

2 

1.0131e  +  00 

0.77 

0.71 

0.77 

0.10 

4 

2.3267e  —  01 

0.76 

0.97 

0.76 

0.73 

8 

5.6522e  —  02 

0.75 

0.99 

0.75 

0.93 

16 

1.4022e  —  02 

0.75 

1.00 

0.75 

0.97 

32 

3.4986e  -  03 

0.75 

1.00 

0.75 

0.98 

64 

8.7422e  -  04 

0.75 

1.00 

0.75 

0.99 

128 

2.1853e  —  04 

0.75 

1.00 

0.75 

1.00 

Conv.  Order  =  2.00 


9i  =  66.987,  g2  =  -2.2807, 


Ne 

|#1  Ql,h,p\ 

V2h,p 

r]h,p+ 1 

ri2h,p 

Vh,p+ 1 

2 

1.1787e  —  01 

0.66 

0.48 

0.66 

2.99 

4 

1.0964e  -  02 

0.15 

1.53 

0.15 

1.00 

8 

6.3700e  -  04 

0.18 

1.98 

0.18 

1.17 

16 

3.8167e  -  05 

0.21 

2.14 

0.21 

3.37 

32 

2.3581e  -  06 

0.22 

2.19 

0.22 

7.82 

64 

1.4696e  —  07 

0.22 

2.22 

0.22 

16.7 

128 

9.1781e  —  09 

0.22 

2.23 

0.22 

34.4 

Conv.  Order  =  4.00 


Ne 

l?2  -  <?2,h,p| 

V2h,p 

Vh,p+ 1 

T]2h,p 

Vh,p+ 1 

2 

5.1806e  -  03 

0.94 

0.95 

0.94 

0.86 

4 

2.5594e  -  04 

0.94 

1.11 

0.94 

0.98 

8 

1.3941e  —  05 

0.94 

1.07 

0.94 

0.95 

16 

8.4825e  -  07 

0.94 

1.02 

0.94 

0.84 

32 

5.2712e  —  08 

0.94 

1.01 

0.94 

0.66 

64 

3.2899e  -  09 

0.94 

1.00 

0.94 

0.31 

128 

2.0555e  —  10 

0.94 

1.00 

0.94 

0.37 

Conv.  Order  =  4.00 


Ne 

\<ls  -  93, h, p| 

Ij2h,p 

nh,p+ 1 

f]2h,p 

Vh,p+1 

2 

1.6235e  —  01 

0.99 

0.79 

0.99 

0.99 

4 

1.1652e  —  02 

0.95 

0.90 

0.95 

1.04 

8 

6.7858e  -  04 

0.94 

0.98 

0.94 

1.10 

16 

4.1144e  —  05 

0.94 

1.00 

0.94 

1.21 

32 

2.5510e  —  06 

0.94 

1.00 

0.94 

1.43 

64 

1.5912e  —  07 

0.94 

1.00 

0.94 

1.87 

128 

9.9398e  -  09 

0.94 

1.00 

0.94 

2.74 

Conv.  Order  =  4.00 


Ne 

\QA  Q4,h,p\ 

r)2h,p 

Vh^+l 

H2h,p 

Vh,p+ 1 

2 

6.5408e  -  03 

0.55 

2.51 

0.55 

38.1 

4 

3.1507e  -  03 

0.93 

0.90 

0.93 

3.05 

8 

2.2060e  -  04 

0.94 

0.97 

0.94 

3.05 

16 

1.4141e  —  05 

0.94 

0.99 

0.94 

6.90 

32 

8.8924e  -  07 

0.94 

1.00 

0.94 

14.8 

64 

5.5663e  -  08 

0.94 

1.00 

0.94 

30.5 

128 

3.4803e  -  09 

0.94 

1.00 

0.94 

62.1 

Conv.  Order  =  4.00 


93  =  149.51,  94  =  -66.795 


Table  6:  Test  2.  Errors  and  effectivity  indexes  772/1, p ,  Vh,p+i ,  V2h,p  and  r]h,p+i  for  the 
output  functionals  q±,  (72,  <73  and  <74  with  Ne  £  {2r}J_2  and  p  =  1,2;  the  asymptotic 
convergence  orders  and  “exact”  values  of  the  output  functionals  are  reported. 
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P  =  1 


p  =  2 


Ne 

Iqi  Ql,h,p\ 

fj2h,v 

Vh,p+ 1 

q2h,p 

Vh,  p+1 

2 

1.0620e  +  00 

0.11 

0.02 

0.11 

0.01 

4 

1.8859e  —  01 

1.10 

0.85 

1.10 

0.75 

8 

4.6160e  —  03 

9.65 

2.33 

9.65 

6.96 

16 

6.6637e  -  03 

0.62 

2.25 

0.62 

0.85 

32 

1.7265e  —  03 

0.58 

2.62 

0.58 

1.07 

64 

4.3536e  —  04 

0.58 

2.71 

0.58 

1.12 

128 

1.0908e  —  04 

0.57 

2.74 

0.57 

1.14 

Conv.  Order  =  2.00 


Ne 

\l2  -  ?2,h,j,| 

7/2  h,p 

Wh,p+1 

V2h,p 

Vh,p+1 

2 

3.6770e  +  00 

0.24 

0.03 

0.24 

0.02 

4 

1.2626e  —  01 

3.00 

0.46 

3.00 

0.57 

8 

6.5856e  —  02 

1.03 

1.42 

1.03 

1.47 

16 

4.4395e  -  03 

0.78 

1.51 

0.78 

1.84 

32 

1.0978e  —  03 

0.76 

1.11 

0.76 

1.61 

64 

2.7624e  —  04 

0.75 

1.03 

0.75 

1.56 

128 

6.9183e  —  05 

0.75 

1.01 

0.75 

1.55 

Conv.  Order  =  2.00 


Ne 

1©  -  Q3,h,p\ 

rl2h,p 

r)h,p+ 1 

q2h,p 

Vh,p+1 

2 

8.3879e  -  01 

0.62 

0.09 

0.62 

0.06 

4 

5.9996e  —  02 

0.15 

0.96 

0.15 

0.27 

8 

2.7870e  —  02 

0.79 

1.13 

0.79 

0.91 

16 

6.4018e  —  03 

0.75 

1.06 

0.75 

0.80 

32 

1.6060e  —  03 

0.75 

1.02 

0.75 

0.76 

64 

4.0219e  —  04 

0.75 

1.00 

0.75 

0.75 

128 

1.0059e  —  04 

0.75 

1.00 

0.75 

0.75 

Conv.  Order  =  2.00 


Ne 

\Q4  Q4,h,p\ 

V2h,p 

Vh,p+ 1 

V2h,p 

Vh,p+1 

2 

1.4515e  +  00 

0.03 

0.01 

0.03 

0.01 

4 

3.1567e  —  01 

1.29 

0.63 

1.29 

0.53 

8 

3.4253e  —  02 

0.34 

0.44 

0.34 

0.13 

16 

1.7935e  —  02 

0.74 

1.01 

0.74 

0.83 

32 

4.5925e  —  03 

0.75 

1.01 

0.75 

0.82 

64 

1.1556e  —  03 

0.75 

1.00 

0.75 

0.81 

128 

2.8938e  —  04 

0.75 

1.00 

0.75 

0.81 

Conv.  Order  =  2.00 


q\  =  1.1153,  q2  =  -2.7978, 


Ne 

|<?1  Ql,h,p\ 

V2h,p 

Vh,p+1 

V2h,p 

Vh,p+ 1 

2 

7.1760e  -  02 

0.50 

1.24 

0.50 

1.26 

4 

4.5158e  —  02 

3.33 

2.07 

3.33 

2.37 

8 

3.2464e  -  02 

0.75 

0.86 

0.75 

0.87 

16 

2.4634e  -  04 

0.24 

2.41 

0.24 

1.36 

32 

1.3965e  -  05 

0.03 

2.05 

0.03 

1.05 

64 

8.1306e  -  07 

0.01 

2.00 

0.01 

0.82 

128 

4.9573e  -  08 

0.01 

1.99 

0.01 

0.41 

Conv.  Order  =  4.04 


Ne 

1 92  -  <?2,7i,p| 

T]2h,p 

r)h,p+ 1 

V2h,p 

Vh,p+ 1 

2 

1.3917e—  01 

0.23 

0.66 

0.24 

0.92 

4 

1.0903e  -  01 

1.83 

0.69 

1.83 

0.84 

8 

4.2557e  -  02 

1.15 

0.61 

1.15 

0.97 

16 

3.2737e  -  04 

0.98 

1.37 

0.98 

1.13 

32 

8.3927e  -  06 

0.95 

1.26 

0.95 

1.10 

64 

3.8305e  -  07 

0.94 

1.08 

0.94 

0.08 

128 

2.2432e  -  08 

0.94 

1.02 

0.94 

2.00 

Conv.  Order  =  4.10 


Ne 

1*73  -  ?3,/i,pl 

H2h,p 

Vh,p+ 1 

H2h,p 

Vh,p+ 1 

2 

4.1045e  —  02 

0.28 

0.49 

0.28 

0.26 

4 

3.0216e  —  02 

1.13 

0.67 

1.13 

0.87 

8 

1.4802e  -  03 

1.50 

0.17 

1.50 

0.70 

16 

8.7423e  -  05 

0.96 

1.24 

0.96 

0.69 

32 

3.8437e  -  06 

0.94 

1.08 

0.94 

1.00 

64 

2.1475e  —  07 

0.94 

1.02 

0.94 

2.13 

128 

1.3055e  —  08 

0.94 

1.01 

0.94 

4.44 

Conv.  Order  =  4.04 


Ne 

\Q4  Q4,h,p\ 

n2h,p 

Vh,p+ 1 

n2h,p 

Vh,p+ 1 

2 

2.1004e  -  02 

1.00 

3.08 

0.98 

5.55 

4 

1.3415e  —  03 

58.2 

44.3 

58.2 

46.4 

8 

4.1508e  —  02 

1.04 

0.96 

1.04 

1.10 

16 

5.5119e  —  04 

0.95 

1.20 

0.95 

1.14 

32 

2.8052e  -  05 

0.94 

1.04 

0.94 

1.05 

64 

1.6060e  -  06 

0.94 

1.01 

0.94 

1.35 

128 

9.7772e  -  08 

0.94 

1.00 

0.94 

2.02 

Conv.  Order  =  4.04 


g3  =  -1.5349,  94  =  0.63397 


Table  7:  Test  3.1.  Errors  and  effectivity  indexes  772/1, P,  Vh,p+i ,  772/1, P  and  77/1, p+i  for  the 
output  functionals  gi,  g2,  g3  and  g4  with  Ne  £  {2r}J_2  and  p  =  1,2;  the  asymptotic 
convergence  orders  and  “exact”  values  of  the  output  functionals  are  reported. 
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P  =  1 


p  =  2 


Ne 

Iqi  Ql,h,p\ 

V2h,p 

Vh,p+ 1 

92  h,p 

f]h,p+l 

2 

2.4354e  —  01 

0.11 

0.07 

0.11 

0.01 

4 

8.8545e  —  02 

0.18 

0.76 

0.18 

0.60 

8 

3.1482e  —  02 

0.14 

1.30 

0.14 

0.56 

16 

9.9188e  —  03 

0.03 

1.80 

0.03 

0.92 

32 

2.5746e  -  03 

0.02 

1.92 

0.02 

0.98 

64 

6.4948e  -  04 

0.02 

1.95 

0.02 

0.99 

128 

1.6274e  —  04 

0.02 

1.96 

0.02 

1.00 

Conv.  Order  =  2.00 


Ne 

|92  “  92,h,pl 

V2h,p 

Wh,p+1 

V2h,p 

Vh,p+ 1 

2 

1.3614e  —  01 

0.79 

0.81 

0.80 

0.81 

4 

1.2298e  —  02 

5.43 

2.06 

5.43 

0.35 

8 

4.5826e  —  02 

1.10 

1.13 

1.10 

1.42 

16 

1.2723e  —  03 

0.83 

1.65 

0.83 

0.79 

32 

2.6283e  —  04 

0.77 

1.25 

0.77 

0.95 

64 

6.3235e  —  05 

0.75 

1.07 

0.75 

0.96 

128 

1.5679e  —  05 

0.75 

1.02 

0.75 

0.97 

Conv.  Order  =  2.00 


Ne 

\l3  ~  93, A, pi 

V2h,p 

Vh,p+ 1 

V2h,p 

Vh,p+ 1 

2 

5.0522e  —  01 

0.73 

0.37 

0.73 

0.36 

4 

1.6131e  —  01 

0.57 

0.72 

0.57 

0.53 

8 

6.6343e  -  02 

0.69 

1.03 

0.69 

0.93 

16 

1.9524e  —  02 

0.74 

1.05 

0.74 

1.02 

32 

5.0777e  -  03 

0.75 

1.02 

0.75 

1.00 

64 

1.2814e  —  03 

0.75 

1.00 

0.75 

1.00 

128 

3.2110e  —  04 

0.75 

1.00 

0.75 

1.00 

Conv.  Order  =  2.00 


Ne 

\Q4  Q4,h,p\ 

V2h,p 

Vh,p+ 1 

V2 h,p 

Vh,p+1 

2 

7.3282e  —  02 

0.98 

0.36 

0.98 

0.38 

4 

5.6749e  -  02 

1.09 

1.04 

1.09 

1.11 

8 

1.1139e  —  03 

0.99 

0.18 

0.99 

2.35 

16 

1.8407e  —  03 

0.77 

1.03 

0.77 

0.89 

32 

4.4339e  -  04 

0.76 

1.01 

0.76 

0.95 

64 

1.0988e  —  04 

0.75 

1.00 

0.75 

0.97 

128 

2.7413e  —  05 

0.75 

1.00 

0.75 

0.99 

Conv.  Order  =  2.00 


gi  =  0.32890,  ®  =  0.87998, 


Ne 

|#1  Ql,h,p\ 

V2h,p 

Vh,p+1 

92A,p 

Vh,p+ 1 

2 

1.6649e  -  01 

0.29 

0.21 

0.29 

0.22 

4 

4.7141e  —  02 

0.55 

0.33 

0.55 

0.34 

8 

1.0316e  —  02 

0.24 

0.48 

0.24 

0.32 

16 

3.8222e  -  04 

0.01 

2.37 

0.01 

1.13 

32 

1.7474e  -  05 

0.01 

2.13 

0.01 

0.95 

64 

9.9520e  —  07 

0.01 

2.04 

0.01 

0.79 

128 

6.0669e  -  08 

0.01 

2.01 

0.01 

0.54 

Conv.  Order  =  4.04 


Ne 

1 92  -  92, A, p| 

V2h,p 

T]h,p+1 

V2h,p 

Vh,p+ 1 

2 

8.9991e  —  03 

7.41 

1.84 

7.35 

0.51 

4 

5.0555e  —  02 

2.29 

1.30 

2.29 

1.41 

8 

4.4598e  -  02 

1.04 

1.12 

1.04 

0.96 

16 

3.0119e  —  05 

0.86 

5.93 

0.86 

4.96 

32 

3.2597e  -  06 

0.95 

1.49 

0.95 

4.41 

64 

1.6636e  -  07 

0.94 

1.14 

0.94 

8.74 

128 

9.3557e  -  09 

0.94 

1.04 

0.94 

17.8 

Conv.  Order  =  4.10 


Ne 

193  -  93, A, pi 

H2h,p 

Vh,p+ 1 

H2h,p 

Vh,p+ 1 

2 

3.0211e  —  01 

0.62 

0.61 

0.62 

0.60 

4 

1.1336e  —  01 

0.82 

0.35 

0.82 

0.25 

8 

1.7674e  -  02 

0.93 

0.91 

0.93 

0.97 

16 

7.6475e  -  04 

0.95 

1.26 

0.95 

1.48 

32 

3.4985e  -  05 

0.94 

1.07 

0.94 

1.64 

64 

1.9942e  -  06 

0.94 

1.02 

0.94 

2.21 

128 

1.2161e  —  07 

0.94 

1.00 

0.94 

3.42 

Conv.  Order  =  4.04 


Ne 

\Q4:  Q4,h,p\ 

n2h,p 

Vh,p+1 

H2h,p 

Vh,p+ 1 

2 

7.0421e  -  02 

1.69 

1.08 

1.68 

1.17 

4 

3.2583e  -  02 

1.31 

1.23 

1.31 

1.21 

8 

5.6490e  -  03 

1.05 

1.17 

1.05 

1.10 

16 

3.5804e  -  05 

0.96 

0.98 

0.96 

0.49 

32 

1.4246e  -  06 

0.94 

1.00 

0.94 

3.90 

64 

8.1164e  —  08 

0.94 

1.00 

0.94 

9.46 

128 

4.9944e  -  09 

0.94 

1.00 

0.94 

19.9 

Conv.  Order  =  4.04 


93  =  -0.75134,  94  =  -0.034175 


Table  8:  Test  3.2.  Errors  and  effectivity  indexes  772/1, p,  Vh,p+i,  V2h,p  and  rjh,p+i  for  the 
output  functionals  gi,  g2,  g3  and  g4  with  Ne  £  {2r}J_2  and  p  =  1,2;  the  asymptotic 
convergence  orders  and  “exact”  values  of  the  output  functionals  are  reported. 
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P  =  1 


p  =  2 


Ne 

1^1  Ql,h,p\ 

^2h,p 

Vh,p+1 

7/2/1, p 

Vh,p+1 

2 

2.1779e  +  01 

0.51 

0.68 

0.60 

0.85 

4 

9.9297e  +  00 

0.03 

0.81 

0.02 

0.41 

8 

1.6917e  +  00 

0.35 

2.18 

0.35 

0.85 

16 

4.0822e  —  01 

0.44 

2.54 

0.44 

0.98 

32 

1.0397e  —  01 

0.44 

2.57 

0.44 

0.99 

64 

2.6120e  —  02 

0.44 

2.57 

0.44 

0.99 

128 

6.5381e  —  03 

0.44 

2.58 

0.44 

1.00 

Conv.  Order  =  2.00 


Ne 

\l2  - 

V2h,p 

Wh,p+1 

V2h,p 

Vh,p+ 1 

2 

5.1712e  —  01 

0.96 

1.28 

1.33 

1.04 

4 

8.9262e  —  03 

0.72 

1.31 

0.41 

0.39 

8 

1.3851e  —  03 

0.65 

1.22 

212 

482 

16 

5.9992e  —  04 

0.73 

1.03 

59.0 

230 

32 

1.6827e  —  04 

0.75 

1.01 

4.17 

15.1 

64 

4.3218e  —  05 

0.75 

1.00 

1.18 

5.25 

128 

1.0876e  —  05 

0.75 

1.00 

0.16 

2.01 

Conv.  Order  =  2.00 


Ne 

I©  -  ?3,h,p| 

»/2  h,p 

Vh,p+ 1 

7/2  h,p 

r)h,p+ 1 

2 

1.3144e  +  00 

0.90 

1.05 

1.07 

1.31 

4 

4.3784e  —  01 

0.88 

0.96 

0.88 

1.11 

8 

8.5190e  —  02 

0.76 

1.06 

0.77 

1.08 

16 

2.1252e  —  02 

0.75 

1.03 

0.76 

1.04 

32 

5.3832e  —  03 

0.75 

1.01 

0.76 

1.01 

64 

1.3504e  —  03 

0.75 

1.00 

0.75 

1.00 

128 

3.3788e  —  04 

0.75 

1.00 

0.75 

1.00 

Conv.  Order  =  2.00 


Ne 

\<Il  ~  Ql,h,p  | 

^]2h,p 

r)h,p+ 1 

^]2h,p 

Vh,p+ 1 

2 

7.6587e  -f  00 

0.37 

0.42 

0.42 

0.29 

4 

5.8921e  —  02 

2.74 

12.0 

3.90 

11.1 

8 

8.9579e  -  02 

0.14 

3.70 

0.14 

1.81 

16 

9.3778e  -  03 

0.34 

1.88 

0.34 

1.16 

32 

4.6134e  -  04 

0.27 

1.77 

0.27 

1.03 

64 

2.7059e  -  05 

0.25 

1.74 

0.25 

1.01 

128 

1.6646e  -  06 

0.25 

1.74 

0.25 

1.00 

Conv.  Order  =  4.02 


Ne 

l<72  -  ©,h,p| 

?\2h,p 

Vh,p+1 

T]2h,p 

Wh,p+1 

2 

1.3746e  -  01 

0.95 

0.46 

3.85 

6.50 

4 

8.9098e  -  03 

0.93 

0.86 

0.14 

2.33 

8 

6.6169e  -  04 

0.94 

1.11 

9.73 

18.9 

16 

4.2862e  -  05 

0.94 

1.03 

2.73 

8.03 

32 

2.5470e  -  06 

0.94 

1.01 

0.55 

2.65 

64 

1.5723e  -  07 

0.94 

1.00 

0.21 

0.94 

128 

9.7966e  -  09 

0.94 

1.00 

0.57 

0.13 

Conv.  Order  =  4.00 


Ne 

1©  -  ©,h,p| 

fj2h,p 

Vh,p+ 1 

fj2  h,p 

Vh,p+ 1 

2 

2.7769e  -  01 

1.08 

0.19 

1.22 

0.32 

4 

6.1107e  —  03 

0.42 

2.50 

0.34 

0.86 

8 

3.8213e  —  03 

0.93 

1.87 

0.94 

1.64 

16 

3.0891e  -  04 

0.95 

1.15 

0.95 

1.13 

32 

1.5794e  -  05 

0.94 

1.04 

0.94 

1.02 

64 

9.3709e  -  07 

0.94 

1.01 

0.94 

0.99 

128 

5.7814e  —  08 

0.94 

1.00 

0.94 

0.99 

Conv.  Order  =  4.02 


©  =  45.293,  q2  =  1.5882, 


©  =  -3.9254 


Table  9:  Test  4.  Errors  and  effectivity  indexes  772/1, p,  Vh,P+i,  772/1, p  and  r)h,P+i  for  the 
output  functionals  <71,  qi  and  <73  with  Ne  6  {2r’}J=2  and  p  =  1,2;  the  asymptotic 
convergence  orders  and  “exact”  values  of  the  output  functionals  are  reported. 


5.2  Discussion 

Based  on  the  results  reported  in  Tables  3-9,  we  extrapolate  the  following  considerations 
for  the  error  estimators  A t^h,P,  ^i,h,P+i,  \,2h,p  and  Aij?l,p+i;  see  Eqs.(31),  (32),  (35) 
and  (36). 

1.  As  expected,  all  the  error  estimators  could  highlight  poor  performances  when 
a  small  number  of  “elements”  Ne  =2,4  is  considered  (see  e.g.  the  effectivity 
indexes  for  Test  1.3  and  gi);  also,  the  effectivity  indexes  largely  changes  from  an 
error  estimator  to  an  other.  Even  if  in  many  cases  the  effectivity  indexes  can 
be  considered  “sufficiently”  close  to  1  (in  particular  for  the  estimators  A %,2h,p> 
Aj^p+i  and  A *,2/j,p),  the  reliability  of  the  estimators  largely  depends  on  the  Test 
problem  for  Ne  ’“small”  ( Ne  =  2,4). 

2.  As  expected,  the  error  estimators  do  not  represent  rigorous  error  bounds  for  the 
errors  on  the  output  functionals  qi~qA  for  all  the  Test  problems. 

3.  For  the  error  estimator  A,  2/i.,p  (Eq.(31)): 
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•  asymptotic  effectivity  indexes  of  values  0.75  and  0.94  are  achieved  with  p  =  1 
and  p  =  2,  respectively,  for  the  output  functionals  (?2~94  (linear  and  nonlinear 
with  the  trigonometric  functions  sine  and  cosine)  in  all  the  tests  considered; 

•  asymptotic  effectivity  indexes  in  the  range  0.20-2.55  are  achieved  for  the 
quadratic  strain  energy  functional  q\  with  rjith,p+i  <  1  in  all  the  Test  prob¬ 
lems  except  for  Test  2  (p  =  1);  for  the  Tests  3.1  [p  =  2)  and  3.2  the  error  on 
qi  is  largely  underestimated  even  asymptotically. 

4.  For  the  error  estimator  Ajylp+1  (Eq.(32)): 

•  the  asymptotic  effectivity  index  =  1.00  is  achieved  with  both  p  =  1 

and  p  =  2  for  the  output  functionals  <72-94  in  all  the  Tests  problems; 

•  asymptotic  effectivity  indexes  bigger  than  one  >  1.00)  are  obtained 

for  the  quadratic  strain  energy  functional  q\  in  all  the  Test  problems;  asymp¬ 
totically,  Ai p^p+i  is  in  the  range  1.33-2.58  for  all  the  Test  problems  and 
differs  for  p  =  1  and  p  =  2. 

5.  For  the  error  estimator  A i)2h,p  (Eq.(35)): 

•  the  behavior  of  this  estimator  is  very  close  to  the  one  of  Aj t2h,p  for  p  =  1  and 
p  =  2  in  many  Test  problems,  except  when  a  small  number  of  “elements”  is 
considered  ( Ne  =  2,4)  or  as  e.g.  in  Test  4  for  q2. 

6.  For  the  error  estimator  Aj^p+i  (Eq.(36)): 

•  asymptotic  effectivity  indexes  are  achieved  only  for  p  =  1  in  a  range  r/ip^+i  = 
0.20-2.18;  these  values  considerably  vary  with  the  Test  problem  and  the 
output  functional  under  consideration  (even  if  linear  <72-94); 

•  the  behavior  of  the  error  estimator  is  quite  unpredictable  for  p  =  2  for  which 
no  asymptotic  effectivity  indexes  are  achieved;  very  large  overestimations  of 
the  errors  (r]i,h,p+i  1)  are  obtained  e.g.  for  Tests  1.1  1.3. 

Based  on  the  previous  considerations  and  by  recalling  Remark  4.3  regarding  the  compu¬ 
tational  costs,  we  conclude  that  the  estimator  Aj^p+i  (31),  which  employs  the  higher 
order  enhanced  dual  solution,  represents  the  most  reliable  and  preferable  error  estima¬ 
tor  among  the  ones  considered  for  the  evaluation  of  the  errors  resulting  from  the  linear, 
quadratic  and  nonlinear  (with  trigonometric  functions  sine  and  cosine)  outputs  and  for 

P  =  1,2. 

If  the  classic  Finite  Element  method  with  linear  Lagrangian  basis  is  used  (see  Re¬ 
marks  4.1  and  4.2),  the  error  estimator  Aj t2h,i  is  preferable  to  A due  to  the  con¬ 
siderations  in  Remark  4.3  regarding  the  computational  costs. 

6  Conclusions 

We  have  proposed  and  studied  goal-oriented  a  posteriori  error  estimators  for  the  evalu¬ 
ation  of  the  errors  on  quantities  of  interest  associated  with  the  solution  of  geometrically 
nonlinear  curved  elastic  rods  with  arbitrarily  large  planar  deflections  under  the  so-called 
Elastica  theory;  for  the  numerical  solution  of  these  problems,  a  Galerkin  formulation 
with  B-spline  basis  of  order  one  and  two  has  been  considered.  We  have  introduced 
goal-oriented  error  estimators  with  higher  order  “enhanced”  primal  and  dual  solutions, 
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in  which  the  enrichment  of  the  original  B-spline  basis  by  means  of  “pure”  fc-refinement 
has  been  used.  These  estimators  have  been  compared  in  several  test  cases  with  the 
most  traditional  ones  based  on  the  mesh  refined  “enhanced”  primal  and  dual  solutions. 
The  numerical  tests  reveal  a  better  performance,  in  terms  of  effectivity  indexes,  for 
the  proposed  error  estimator  based  on  the  higher  order  “enhanced”  dual  solution  with 
respect  to  one  based  on  the  mesh  refined  “enhanced”  dual  solution;  while  asymptotic  ef¬ 
fectivity  indexes  larger  than  one  have  been  obtained  for  the  quadratic  functional  (strain 
energy),  unitary  asymptotic  effectivity  indexes  have  been  achieved  for  the  linear  and 
the  nonlinear  output  functionals  (rotation  and  displacements). 

The  error  estimators  can  eventually  be  adopted  for  the  analysis  of  other  structural 
problems,  such  as  e.g.  the  so-called  Kirchhoff  and  Reissner  Sinro  spatial  beam  theories, 
or  other  more  general  PDEs. 
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